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1 Introduction 

Twenty years ago Car and Parrinello introduced an efficient method to perform 
Molecular Dynamics simulation for classical nuclei with forces computed on the 
"fly" by a Density Functional Theory (DFT) based electronic calculation[l]. Because 
the method allowed study of the statistical mechanics of classical nuclei with many- 
body electronic interactions, it opened the way for the use of simulation methods for 
realistic systems with an accuracy well beyond the limits of available effective force 
fields. In the last twenty years, the number of applications of the Car-Parrinello ab- 
initio molecular dynamics has ranged from simple covalent bonded solids, to high 
pressure physics, material science and biological systems. There have also been ex- 
tensions of the original algorithm to simulate systems at constant temperature and 
constant pressure[2], finite temperature effects for the electrons [3], and quantum 
nuclei [4]. 

DFT is, in principle, an exact theory but the energy functional are treated ap- 
proximately at the level of a self consistent mean field theory for practical purposes. 
Despite recent progress, DFT suffers from well-known limitations, for example, ex- 
cited state properties such as optical gaps and spectra are generally unreliable. DFT 
shows serious deficiencies in describing van der Waals interactions, non-equilibrium 
geometries such as reaction barriers, systems with transition metals and/or cluster 
isomers with competing bonding patterns [5, 6]. As a consequence, current ab-initio 
predictions of metallization transitions at high pressures, or even the prediction of 
phase transitions are often only qualitative. Hydrogen is an extreme case [7, 8, 9] but 
even in silicon the diamond//3-tin transition pressure and the melting temperature are 
seriously underestimated [10]. 

An alternative route to the ground state properties of a system of many electrons 
in presence of nuclei is the Quantum Monte Carlo method (QMC) [11, 6]. QMC 
methods for bosons are "exact " meaning that all systematic eiTors are under con- 
trol and can be reduced as much as desired with a computational cost growing as a 
power of the number of particles. However for fermions the "sign problem" makes 
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a direct extension of QMC unstable and one has to resort to the "fixed node approx- 
imation" for practical calculations[l 1, 6|. Over the years, the level of accuracy of 
the fixed node approximation for simple homogeneous systems, such as '^ife and the 
electron gas, has been systematically improved [12, 13, 14]. In more complex, inho- 
mogeneous situations such as atoms, molecules and extended systems of electrons 
and nuclei, progress have been somewhat slower. Nonetheless, in most cases, fixed- 
node QMC methods have proved to be more accurate than mean field methods (HF 
and DFT)[6]. Computing ionic forces with QMC to replace the DPT forces in the 
ab-initio MD, is more difficult and a general and efficient algorithm is still missing. 
Moreover, the computer time required for a QMC estimate of the electronic energy 
is, in general, more than for a corresponding DFT-LDA calculation. These problems 
have seriously limited the development of an ab-initio simulation method based on 
the QMC solution of the electronic problem "on the fly". 

In recent years, we have developed a different strategy based entirely on the 
Monte Carlo method, both for solving the electronic problem, and for sampling the 
ionic configuration space[15, 16]. The new method, called the Coupled Electron-Ion 
Monte Carlo method (CEIMC) relies on the Born-Oppenheimer approximation for 
treating finite temperature ions coupled with ground state electrons. A Metropolis 
Monte Carlo simulation of the ionic degrees of freedom (represented either by clas- 
sical point particles or by path integrals) at fixed temperature is performed based on 
the electronic energies computed during independent ground state Quantum Monte 
Carlo calculations. CEIMC has been appUed, so far, to high pressure metalUc hydro- 
gen where it has found quite different effects of temperature than CPMD with Local 
Density Approximation (LDA) forces[17]. In these lecture notes, we present the the- 
oretical basis of CEIMC. We start by describing, in some detail, the ground state 
QMC methods implemented in CEIMC, namely the Variational Monte Carlo and the 
Reptation Quantum Monte Carlo methods. We then describe the fixed node or re- 
stricted paths approximation necessary to treat fermions and the fixed phase method 
used to perform the average over boundary conditions needed for metallic system. In 
the subsequent section we describe how to implement the Metropolis algorithm when 
the energy difference has a statistical noise and discuss efficient strategies for energy 
differences within the CEIMC. The next section is devoted to describe the method 
to treat quantum mechanical protons and to integrate efficiently this new difficulty 
within CEIMC. Finally, we briefly review some CEIMC results for high pressure hy- 
drogen and compare with existing CPMD results. Conclusions and perspectives for 
future developments are collected in the last section. 

2 The electronic ground state problem 

Let us consider a system of Np nuclei and A^e electrons in a volume V described by 
the non relativistic hamiltonian 




(1) 



THE COUPLED ELECTRON-ION MONTE CARLO METHOD 



3 



where -/V = N^ + Np, and 2;,, mj, fj represent the charge, mass and position operator 
of particle i respectively and — fi^ /2mi. Let us denote with R = (ri, • • • , rjVe) 
and S = (rjVe+ij • • • ) ''iv) the set of coordinates of all electrons and nuclei respec- 
tively. We restrict the discussion to unpolarized systems, i. e. systems with a van- 
ishing projection of the total spin along a given direction, say Sz ~ 0. Since the 
Hamiltonian does not flip spins, we can label electrons from 1 to Ne/2 as up spin (t) 
and electrons from Nf./2 + IX.0 as down spin (|). 

Within the Born-Oppenheimer approximation, the energy of the system for a 
given nuclear state S is the expectation value of the hamiltonian H over the corre- 
sponding exact ground state \^o{S) > 

Ebo{S) = < MS)\H\MS) > (2) 

which is a SNg dimensional integral over the electronic coordinates in configura- 
tional space 

EBoiS)^ J dR<P*o{R\S)H{R,S)MR\S)= J dR\<Po{R\S)\^ El{R\S) (3) 
with the local energy defined as 

Since (^q{R\S) is normalized and |<?o(^|'S')|^ > everywhere, the 3A^e dimensional 
integral in eq. (3) can be performed by standard Metropolis Monte Carlo by generat- 
ing a Markov process which sample asymptotically |^o(^|5')|^. Expectation values 
of any observable can be computed along the same Markov chain. In this respect, 
computing the properties of a many-body quantum system is similar to performing a 
MC calculation for a classical system. The square modulus of the ground state wave 
function plays the role of the classical Boltzmann distribution. An important quantity 
in what follows is the measure of the energy fluctuations for a given wave function. 
This can be defined by the variance of the local energy 

dR(^HMR\S)y -EloiS) (5) 



Note that for any exact eigenfunction of the hamiltonian, the local energy El{R\S) 
does not depend on the electronic configuration R and is equal to the correspond- 
ing eigenvalue. This implies that the variance vanishes. This is known as the zero 
variance principle in Quantum Monte Carlo. 



2.1 Variational QMC 



The problem is to get a computable expression for the ground state wave function 
without solving the Schrodinger equation for the many body hamiltonian of eq. (1), 
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obviously an impossible task for any non trivial system. As usual in many body 

problems, we can resort to the variational principle which states that the energy of 
any proper trial state \^t{S) > will be greater or equal to the ground state energy 

p IdR1'UR\S)H{R,S)W\S) 
EsoiS) < EviS) = jdji^.^ji\sPr{R\S) 

A proper trial wave function must satisfy the following requirements 

• it has to have the right symmetry under particle permutation: ]i'T{PR\S) = 
(— 1)^^t{R\S), where P is the permutation operator for electrons of same spin. 

• the quantity H'Ft needs to be well defined everywhere which implies that both 
iZ't and Vlf'T must be continuous whenever the potential is finite, including at 
the periodic boundaries. 

• the integrals / (ii?|!Z'Tp and / dR^^H^r must exist. Furthermore for a Monte 
Carlo evaluation of the variance the integral / dR{H\l/T)'^is also required to 
exist. 

For a given trial function, it is essential to show analytically that these properties 
hold everywhere, in particular at the edge of the periodic box and when two particles 
approach each other (where generally the potential diverges). Otherwise, either the 
upper bound property is not guaranteed or the Monte Carlo error estimates are not 
valid. 

The strategy in Variational Monte Carlo (VMC) is therefore to pick a proper form 
for a trial wave function based on physical insight for the particular system under 
study. In general, a number of parameters {ai , • • • , a^.) will appear in the wave func- 
tion to be treated as variational parameters. For any given set of {a} the Metropolis 
algorithm is used to sample the distribution 

and the electronic properties are then computed as averages over the generated 
Markov chain 

E^{S,{a})=<EL{R\S,{a})> (8) 

where the superscript T as been explicitly written to remember that the variational 

energy depends in general on the chosen analytical form and, for any given form, on 
the numerical values of the variational parameters {a}. 

Because of the zero variance principle stated above, the fluctuations in the local en- 
ergy are entirely due to inaccuracies of the trial function for the particular config- 
urations generated during the MC run. As the trial wave function approaches the 
exact eigenfunction (everywhere in configuration space!) the fluctuations decrease 
and the variational estimate of the energy converges more rapidly with the number 
of MC steps. At the same time, the estimate converges to the exact energy. This is at 
variance with classical Monte Carlo where fluctuations are induced by temperature. 

The variational method is very powerful, and intuitively pleasing. One posits a 
form of the trial function and then obtains an upper bound for the energy. In contrast 
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to other theoretical methods, no further approximations are made. The only restric- 
tion on the trial function is to be computable in a reasonable amount of time. 

One of the problems with VMC is that it favors simple states over more com- 
plicated states. As an example, consider the liquid-solid transition in helium at zero 
temperature. The solid wave function is simpler than the liquid wave function be- 
cause in the solid the particles are localized so that the phase space that the atoms 
explore is much reduced. This biases the difference between the Uquid and solid 
variational energies for the same type of trial function, ( e.g. a pair product form, see 
below) since the solid energy will be closer to the exact result than the hquid. Hence, 
the transition density will be systematically lower than the experimental value. An- 
other illustration is the calculation of the polarization energy of liquid ''He. The 
wave function for fully polarized helium is simpler than for unpolarized hehum be- 
cause antisymmetry requirements are higher in the polarized phase so that the spin 
susceptibility computed at the pair product level has the wrong sign! 

The optimization of trial functions for many-body systems is time consuming, 
particularly for complex trial functions. The dimension of the parameter space in- 
creases rapidly with the complexity of the system and the optimization can become 
very cumbersome since it is, in general, a nonlinear optimization problem. Here we 
are not speaking of the computer time, but of the human time to decide which terms 
to add, to program them and their derivatives in the VMC code. This allows an ele- 
ment of human bias into VMC; the VMC optimization is more likely to be stopped 
when the expected result is obtained. The basis set problem is still plaguing quantum 
chemistry even at the SCF level where one only has 1-body orbitals. VMC shares 
this difficulty with basis sets as the problems get more complex. 

Finally, the variational energy is insensitive to long range order. The energy is 
dominated by the local order (nearest neighbor correlation functions). If one is try- 
ing to compare the variational energy of a trial function with and without long range 
order, it is extremely important that both functions have the same short-range flex- 
ibility and both trial functions are equally optimized locally. Only if this is done, 
can one have any hope of saying anything about the long range order. The error in 
the variational energy is second order in the trial function, while any other property 
will be first order. Thus variational energies can be quite accurate while correlation 
functions are not very accurate. 

As a consequence, the results typically reflect what was put into the trial function. 
Consider calculating the momentum distribution. Suppose the trial function has a 
Fermi surface. Then the momentum distribution will exhibits a discontinuity at kf 
signaling the presence of a Fermi surface. This does not imply that the true wave 
function has a sharp Fernoi surface. 

2.2 Reptation Quantum Monte Carlo 

It is possible to go beyond VMC by a number of related methods known as Pro- 
jection Monte Carlo Methods. The general idea is to chose a trial function which 
has a non negligible overlap with the ground state wave function (in general, it is 
enough to require the right symmetry) and to apply a suitable projection operator 
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which zeros out all the components of the trial wave function from the excited states 
in the Hilbert space of the system. We will limit the description to the method imple- 
mented in CEIMC, namely the Reptation Quantum Monte Carlo [18] or Variational 
Path Integral [19, 20]. For the discussion of other projection QMC methods such as 
Diffusion Monte Carlo (DMC) and Green Function Monte Carlo (GFMC), we refer 
to the specialized literature[6, 11]. 

Let us define {<l>i, Ei} as the complete set of eigenfunction and eigenvalues of the 
hamiltonian H in eq.(l). Any trial state can be decomposed in the eigenstate basis: 

\^T>=Y,C^\<Pi> (9) 
i 

where Ci is the overlap of the trial state with the i*^ eigenstate. Let us consider the 
application of the operator e~*^ onto this state 

i 

with the initial state |!f'(0) >= \\1/t >■ Here i is a control parameter with dimension 
of inverse energy and we will call it "time" since it plays the role of imaginary time 
in the Bloch equation (see below). All excited states will be zeroed exponentially fast 
with increasing t, the rate of the convergence to the ground state depending on the 
energy gap between the ground state and the first excited state non-orthogonal to the 
trial function. The total energy as function of time is defined as 

_ < l^it/2)\HMt/2) > _ < ^T\e-i"He-i"\^T > 
^' <^{t/2)\^{t/2)> < !i^T|e-*^|a^T > 

Similar to a thermal partition function, let us define the generating function of the 
moments of if as 

Z{t) WT\e-*"\<I'T > . (12) 
The total energy at time t is simply the derivative of the logarithm of Z{t) 

E{t) = -^nZit) (13) 
and the variance of the energy is the second derivative 

al{t)=<{H-E{t)f>=-^^E{t)>0 (14) 

which is non-negative by definition. This implies that the energy decreases mono- 
tonically with time. The ground state is reached at large time (much larger than the 
inverse gap) and 

lim E{t) = Eo (15) 
limo-2(f)=0 (16) 

t^OO 
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The last relation is the generalization of the zero variance principle in Projection 
Monte Carlo. 

For observables A which do not commute with H, for instance correlation functions, 
the average at "time" t, defined as in eq. (11), takes the following form in configura- 
tional space ^ 

A{t) =< A>t=-^j dRrdR^dRidRi < ^t\Ri > p{RuR2\\) 

<R2\A\R3> p{R3,R4\^)<R4^T> (17) 

where Ri represent the set of all electronic coordinates and p{R, R', t) is the thermal 
density matrix of the system at inverse temperature t 

p{R,R',t)=<R\e-*"\R' > (18) 
Similarly, the expression of Z{t) in configurational space is 

Z{t) = j dRidR2 < ^t\Ri > p{Ri,R2,t) < R2WT > (19) 

Thus, in order to compute any average over the ground state we need to know the 
thermal density matrix at large enough "time". Obviously, its analytic form for any 
non-trivial many-body system is unknown. However, at short time (or high tempera- 
ture) the system approaches its classical limit and we can obtain approximations. Let 
us first decompose the time interval t in M smaller time intervals, r = t/M 

. M-l 

p{R,R',t)=<R\e-^^"^''\R' >= dRi---dRM-i J] < Rk-i\e-^^\Rk > 

■' fe=i 

(20) 

with the boundary conditions: Rq = R and Rm = R' on the paths. For M large 
enough, we can apply the Trotter factorization to get an explicit form for the short 
time propagator. The simplest factorization, known as the "primitive" approximation, 
consists of ignoring the commutator of the kinetic and potential operators 

p{Rk-i,Rk,r) =< Rk-i\e-^"\Rk >^< Rk-i\e-^^\Rk > e"* 

(21) 

A more accurate, but also more complex form, will be discussed later in the section 
on Path Integral Monte Carlo. Note that we have symmetrized the primitive form in 
order to reduce the systematic error of the factorization [19]. The exphcit form of the 
kinetic propagator is the Green's function of the Bloch equation of a system of free 
fermions [21, 19], i.e. a diffusion equation in configurational space 

< i?._i|e-^|ii, >= f^J e ^ (22) 



^ the expression gets slightly easier for observables diagonal in configurational space: < 

R\A\R' >= A{R)d{R- R'). 
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and therefore we get 



p{R,R',t) = YldRk 



M-l 



fe=l 



fe=i (47rAr) 



3N/2 



(«m)1 



(23) 

In the continuous limit (M — > oo, r — > 0, t = Mr =const.) it becomes the 
Feynman-Kac formula [19] 



p{R, R', t) = ( exp 



dTViRir)) 



(24) 



RW 



where < • • • >rw indicate a path average over gaussian random walks R{t) starting 
at -R(O) = R and ending at R{t) = i?' in a time t. 

We have, in principle, developed a scheme for Monte Carlo calculations of ground 
state averages of a general quantum system. 

However, this scheme has a serious problem of efficiency which prevents its use 
for any non-trivial system. At the origin of the problem are the wild variations of the 
potential V{R) in configuration space. There are cases like electron-proton systems 
where the potential is not bounded and therefore the primitive approximation of the 
propagator is not stable for any finite time step r. However, even with well behaved 
effective potentials (like in Helium for instance) the large fluctuations of the poten- 
tial energy would require a very small time step in order to observe convergence of 
the averages to their exact value. Moreover, the efficiency will degrade rapidly with 
more particles. The problem was recognized in the early days of QMC and the rem- 
edy introduced by Kalos in 1974. In the community of Ground State QMC it goes 
under the name of "importance sampling" (IS). A different strategy is applied in the 
PIMC community. We now describe importance sampling, not in the original form 
as introduced by Kalos in Green Function Monte Carlo, but following a recent de- 
velopment by Baroni and Moroni in the framework of the Reptation QMC [18]. 
In VMC, a good trial function should have a local energy almost constant in con- 
figuration space. Let us assume to know such function ]Pt- Let us then rewrite the 
hamiltonian H in terms of a new fictitious hamiltonian H 



H = n + EL{R) 



where 



H = X 



2 . V^tf^j. 



-V^-h 



El{R) = V{R)-\ 



\I/t 
't 



(25) 

(26) 
(27) 



We can now factorize the short time propagator in a different way (confr. eq. (21) ) 

p{Rk-i,Rk,T) Rk_,\e-^^\Rk > e-5[^-(«'=)+^-(«'=-01 (28) 
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In this new form, the widely oscillating potential energy is replaced by the local en- 
ergy which is much smoother for an accurate We need to find the short time 
propagator of the importance sampling hamiltonian which is nothing but the solu- 
tion of the corresponding Bloch equation [21, 19] 

-dtp.s {R. R\t)= Hp,, (i?, R', t) (29) 
p,siR,R\0) = S{R- R') (30) 

It is not difficult to show by direct substitution, that the short time solution of this 
equation is 

^T{Rk-i)( I \^ ( {Rk-Rk-i-2XTFk-iy 
PM-uRk,r) = ^^[^^) exp| — 

(31) 

if we make the short time approximation [l + r (V-F + V^-F) ~ l] . In these ex- 
pressions, the drift force is defined as Ffc = F{Rk) = 2V In ^T{Rk)- 
This form of the short time propagator does not satisfy an important property of den- 
sity matrices, namely the symmetry under exchange of the two legs, R and R'. We 
can remedy by taking the symmetrized density matrix as short time propagator 



p%{Rk-i,Rk,T) = [p,s{Rk-i, Rk,T)p,s{Rk, Rk-i,T)]^ = 



^-Ls{Rk-i,Rk,r) 



3N 

(47rAr) ^ 

(32) 

where the expression for the symmetrized link action is 

T fT? D ^\ {Rk-Rk-lf >^T /p2 , p2 , {Rk- Rk-i) ■ {Fk - -Pfc-l) 
Ls[Hk-i,Hk,T) = — \--^{i'k+^k-i)-^ 2 

(33) 

Using eqs. (32), (33) we obtain the propagator at any time t as 



piR,R',t) = f Y[dRk 

J 1 1 



■.3N/2 

.fe=i (47rAr) 



(34) 



In the continuum limit it is the generalized Feynman-Kac formula 

p{R, R',t) = /exp (- / dTEL{R{T))) ) (35) 

\ V ^0 / / DRW 

where < • • • > drw indicate a path average over drifted random walks starting at 
i?(0) = R and ending at R{t) = i?' in a time t. With this form of the density matrix 
the generating function takes the form 

Z{t)= f dRdR'^T{R){e~^o'^^^'^^^^^^^\ ^t{R') (36) 

\ / DRW 



and the average of a generic observable A becomes 
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^ J \ / DRW 

\ / DRW 

with obvious boundary conditions on the path averages. 

A special word on the calculation of the energy is in order. In the last equality in 
eq. (1 1) the hamiltonian operator in the numerator can be pushed either to the left 
or to the right in such a way to operate directly on the trial state. Remembering the 
definition of the local energy, eq. (4), we can write 

E{t) = J dRdR'EL{R)^T{R)p{R,R',tMR') 
^W}J dRdR'^T{R)p{R,R' ,t)^{R')EL{R') 
= ^ < El{R) + El{R') > (38) 

We use the last equality in order to improve the efficiency of the estimator. When 
computing the variance we push one H operator to the left and the other to the right 
to obtain 

a^{t)= J dRdR'EL{R)^T{R)p{R,R',t)^{R')EL{R')- E^{t) (39) 

Then reaching the ground state with vanishing variance means taking paths long 
enough (t large enough) for the correlation between the two ends to vanish. On the 
other hand the VMC method is obtained for t = as p{R, R', 0) = d{R - R'). 



2.3 Fermions 



Up to now we have tacitly ignored the particle statistics and derived the formalism 
as if the particles were distinguishable. As far as the Hamiltonian does not depend 
explicitly on spin, the formalism remains valid for fermions if we consider states 
completely antisymmetric under particle permutation I Pi?i >= {—)^\Ri >. The 
importance sampling hamiltonian H, defined in the same way, is symmetric under 
particle exchange even for antisymmetric trial states. The only place where we need 
care is in the initial condition of the Bloch equation, eq. (30), which must be replaced 
by a completely antisymmetric delta function 

p,, {R, R', 0) = AS{R -R') = ^ ^(-l)^^(ii - PR') (40) 

p 

Since [H, P] = 0, the imaginary time evolution preserves the symmetry and the 
fermion thermal density matrix takes the form 
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pAR,R',t) = -^^(-i)^p,(i?,Pi?',i) (41) 
p 

where p„ is the density matrix of a system of distinguishable particles derived above 
(see eq. (35)). The fermion density matrix between configurations R and R' at "time" 
t is the sum over permutations of the density matrix of distinguishable particles be- 
tween the initial configurations R and the permutation of the final configuration PR', 
multiplied by the sign of the permutation. Each of those density matrices arises from 
the sum over all paths with given boundary conditions in time as expressed by the 
generalized Feyrmian-Kac formula. We can therefore think of a path in the config- 
urational space of distinguishable particles as an object carrying not only a weight 
(given by the exponential of minus the integral of the local energy along the path) 
but also a sign fixed by its boundary conditions in time. The fermion density matrix 
is the algebraic sum over all those paths. While p,-, {R, R',t) > for any R' and t at 
given R, this property obviously does not apply to p^ {R, R',t). This is at the origin 
of the "fermion sign problem" [26]. Briefly, the sign problem arises from the fact 
that the optimal probability to sample the electronic paths is the absolute value of 
the fermion density matrix which, however, is a bosonic density matrix (symmetric 
under particle permutation). With this sampUng, the sign of the sampled paths will 
be left in the estimator for the averages. The normalization of any average will be 
given by the number of sampled positive paths minus the number of sampled neg- 
ative paths. Since the sampUng is bosonic, i.e. synraietric, these two numbers will 
eventually be equal and the noise on any average will blow up for a long enough 
sampling. The fundamental reason behind this pathology is that the Hilbert space of 
any time independent Hamiltonian (with local interactions) can be divided in the set 
of symmetric states, antisymmetric states and states of mixed symmetry. These sets 
are disjoint, e. g. any symmetric state is orthogonal to any antisymmetric one; in prin- 
ciple we cannot extract information for a fermionic state from a bosonic samphng 
[22]. 

A general solution of the fermion sign problem is still unavailable, although in- 
teresting algorithms have been proposed [11]. A class of methods try to build the 
antisymmetry constraint into the propagator, while other methods try to reformulate 
the problem of sampling in the space of antisymmetric wavefunctions (determinants) 
[23, 24, 25]. All these "fermions" methods are still at an early stage and their appli- 
cation has been limited so far to quite small numbers of fermions. The more robust 
and widely used, although approximate, method is the so-called restricted path or 
fixed node method [6, 26] . 

Within the fixed-node method, we need to consider the nodal surfaces of the 
fermion density matrix. For any given configuration R, these are defined by the im- 
plicit equation {R, R',t) — 0, as the locations R' at which the density matrix at 
time t vanishes. The nodal surfaces of the initial configuration R divide the config- 
urational space of R' in regions of positive Pp and regions of negative Pp . In terms 
of individual paths, the nodal locations are hypersurfaces in configurational space on 
which the sum of contributions of the positive and negative paths to the density ma- 
trix vanishes. Since the fermion density matrix satisfies the usual convolution relation 
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[21] 

{R, R', t) = j dR"p, {R, R", t)p, {R", R',t-T) Vr e [0, t] (42) 

the configurations R" belonging to the nodal surface of the initial point R at the 

arbitrary time t will not contribute to the integral and therefore to the density matrix 
at any future time t. Therefore in constructing the fermion density matrix from R to 
R' at time t as sum over signed paths, we can safely disregard all those paths which 
have reached the nodal surface at any previous time t < t. If we define the reach 
of R at time t, T{R, t), as the set of points that can be reached from i? in a time t 
without having crossed the nodal surfaces at previous times, the argument above can 
be formaUzed in the restricted paths identity'* 

pAR,R',t) = ^,Y.^-r{l VYe-'A (43) 

' p \JY{0)=R,Y{t)=PR' J r(R,t) 

where S\Y] represent the action of the generic path Y (see eq. (34) for its discretized 
form). Let us now consider the generating function Z{t) of eq. (19). Using the re- 
stricted paths identity we obtain 

Z{t) = f dRdR'^T{R)^,Y.^-f ( f VY e-^M ) MR') 

J ^ \JY{0)=R,Y{t)=PR' J Y(R,t) 

= [dRdR"MR)([ We-^[^1) -^y^i-fMP~'R") 

J \JYiO)=R,Yit)=R" JriR,t)^- P 

= f dRdR"^T{R) I / 'DY e-^[^l ) ^t{R") (44) 

J \JY{0)=R,Y{t)=R" 



r{R,t) 



where in the last equality we have used the fact that the trial wave function is anti- 
symmetric under particle permutations: {—)^^t{P~^R) = !Z't(-R)- In this last form 
the integrand is always positive since for each R, the functional integral is restricted 
to paths inside its reach so that ^t{R)'I't{R') > (if the nodes of the trial function 
are correct). Therefore using the restricted path identity we have proven that the gen- 
erating function is a positive function at any time t and can be computed considering 
only positive paths which do not cross the nodal surfaces. 

The restricted paths identity is by no means the solution of the sign problem 
since in order to know the nodal surfaces we have to know the density matrix itself. 
However rephrasing the problem in terms of spacial boundary conditions can lead 
to interesting approximate schemes. The nodal surface of the fermion density matrix 
for a system of N interacting particles is a highly non-trivial function in 6N dimen- 
sions and not much is known about it [27, 28]. The approximate method, known as 
fixed node in ground state QMC [6] and as restricted paths in PIMC [26], consists in 



in ref. [27] an alternative proof based on the Bloch equation is provided. 
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replacing the nodal surfaces of the exact fermion density matrix with the nodal sur- 
faces of some trial density matrix. In ground state QMC, it is customary to restrict 
the class of possible nodal surfaces to time independent nodes and, within this class, 
the most reasonable choice is to assume the nodes of the trial wave function !Pt. In 
practice, this step requires a very minor modification of the algorithm: it is enough 
to ensure that 'l/T{Rk-i)^T{Rk) > for any time interval along the sampled paths. 
In the continuous limit this restriction will enforce the restricted path identity. 

The nodal surfaces of the trial wave function divide the configurational space in 
disconnected regions. In order to perform the configurational integral over R and 
R' in eq. (44) it could appear necessary to sample all nodal regions. However, it 
can be proved that the nodal regions of the ground state of any Hamiltonian with a 
reasonable local potential are all equivalent by symmetry (Tihng theorem) [27]. This 
"tiling" theorem ensures that computing in a single nodal region is equivalent to a 
global calculation. 

A further important property of the Fixed node method is the existence of a vari- 
ational theorem: the FN-RQMC energy is an upper bound of the true ground state 
energy Et{oo) > Eq, and the equality holds if the trial nodes coincide with the 
nodes of the exact ground state [6]. Therefore for fermions, even projection meth- 
ods such as RQMC are variational with respect to the nodal positions; the nodes are 
not optimized by the projection mechanism. The "quaUty" of the nodal location is 
important to obtain accurate results. 

In some cases it is necessary to consider complex trial functions, for instance in 
the presence of a magnetic field [29] or in the twist average method to be discussed 
later[30]. In these cases we have to deal with a trial function of the form 



where (p^ {R), a real function, is the configuration dependent phase of the wave func- 
tion. Obviously in VMC no differences arise in the use of complex trial functions 
other than in the estimators for the averages. For instance, the local energy is modi- 
fied (confr. eq. (27)) and contains an imaginary part 



El{R) = V{R) - A^i^ + A (V^J' - zA (2y^,^P^ + VV.) (46) 



Here, we limit the discussion to systems with time-reversal symmetry, i. e. zero mag- 
netic fields. As \1/t approaches the exact ground state, the local energy approaches a 
real constant equal to the ground state energy while the imaginary part of the local 
energy vanishes. For general complex functions, we can split the time independent 
Schrodinger equation into two coupled equations, one for the modulus |#| and one 
for the phase (p of the wave function 



(45) 



{-AV^ + [V{R) + XiVifif] } |f I = E\<P\ 



(47) 



o Vl^l 
VV + 2V(^ • = 



(48) 
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It is therefore natural to take the real part of eq. (46) as energy estimator and, in ad- 
dition to the variance, to monitor the deviation of the imaginary part from zero as an 
indicator of the quality of the trial wave function. 

How do we have to modify the RQMC to work with complex wave functions? 
In the "fixed phase" approximation [29, 31] one keeps the phase (Pj.{R) fixed to 
some analytic form during the calculation, and solves the imaginary time dependent 
Schrodinger equation corresponding to the stationary problem of eq. (47). Even for 
fermions this is a bosonic problem (since the modulus of the wave function must be 
symmetric under particle exchange) with a modified interaction [V^(-R) + A(V<p)^] . 
We can still perform the IS transformation and the formaUsm remains the same if 
the local energy is defined as the real part of eq. (46). Note that the fixed node con- 
straint for real trial functions can be recast into the fixed phase algorithm if we write 
friR) = 7r[l — 6'( )] so that the phase of the trial function changes by TT across 

the nodes at it should. Since the phase is a step function, its gradient is a 5 function 
and provides an infinite contribution in the action of paths crossing the nodes, i.e. a 
vanishing probabihty to cross the nodes. 

In this section we have not explicitly indicated the dependence on the ionic state 
S. In the BO approximation, ions play the role of external fields for the electronic 
system so that their coordinates appear explicitly in the Hamiltonian, in the trial state, 
in the local energy and in drift force for the IS procedure. 

2.4 IVial wave functions for hydrogen 

In this subsection, we describe some general properties of the trial wave functions 
for electronic systems. We will restrict our discussion to the case of a proton-electron 
system, i. e. hydrogen and refer to the Uterature for more complex systems (heavier 
elements)[l 1, 6]. In particular, we will not discuss the use of pseudopotentials in 
QMC and the related trial wave functions which, however, will be an important issue 
in future extensions of CEIMC. 

The Pair Product Trial Function 

The pair product trial wave function is the simplest extension of the Slater determi- 
nant of single particle orbitals used in mean field treatment of electronic systems (HF 
or DFT). This is also the ubiquitous form for trial functions in VMC 



where E = {cri, . . . , <7n^} is the set of spin variables of the electrons, 9k{r, (j\S) is 
the A:*'' spin orbital for the given nuclear configuration and 9k{ri, <Ji\S) is the Slater 
matrix. The additional term Uij {rij ) is the "pseudopotential" or parr correlation fac- 
tor which introduces explicitly the two body correlations into the many body wave 
function. This term is of bosonic nature, i.e. is symmetric under particle exchange. 




(49) 
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while the antisymmetry is ensured by the determinant. Often the general form of 
both 9k and Uij are derived by some appropriate theory and then used in connection 
with some free variational parameters to be optimized. 

Let us discuss first the appropriate form of the "pseudopotential" for Coulomb sys- 
tems. There are important analytical properties for the "pseudopotential" that can be 
easily derived. Consider bringing two particles together and let us examine the dom- 
inant terms in the local energy. In a good trial action, the singularities in the kinetic 
energy must compensate those in the potential energy. The local energy for the two 
particle system is 

eL {Uj ) = e"'^ hij [e"'''^ ] = Vij (ry ) + Ay [V^u^ - ( Vuy f] (50) 

where Xij = Aj + Aj, hij is the two body hamiltonian with the interaction potential 
Vij and the trial wave function exp[— u^]. Spin symmetry has been disregarded as 
well as the trivial term related to the center of mass motion. Therefore the short dis- 
tance behavior of any good form of the pseudopotential should follow the solution 
of the two body Schrodinger equation. For the Coulomb potential the "cusp" condi- 
tion derives from this constraint. Indeed substituting r in eq.(50) and 
zeroing the coefficient of the dominant power of r for r — > provides 



duij 



dr 



Xij{D-l) 



(51) 



It is also important to reproduce the correct behavior at large distances where a 
description in terms of collective coordinates is appropriate. The long-wavelength 
modes are important for the low energy response properties and are also the slowest 
modes to converge in QMC. It is possible to show that, within the Random Phase 
Approximation (RPA), the local energy is minimized by imposing 

ul^ = -^ + VT+^ (52) 

with flfe = 12rs/k'^ and the electron sphere radius is related to the volume per 
electron by v = 47rrg/3 in atomic units [32]. The obtained form of the "pseudopo- 
tential" is correct at short and long distances but not necessarily in between because 
of the approximation. One can improve slightly the quaUty of the VMC results con- 
sidering the form 

Uij{r)=u^^''^{r)-a^,e-^'/< (54) 

with the variational parameter aij,Wij. The additional term preserves the short and 
long distance behavior of the RPA function. This form of the pair trial function in- 
troduced four variational parameters, namely aee, Wee, ciep, Wep. 

We discuss now the choice of the spin orbitals. The spin-orbitals are conceptually 
more important than the pseudopotential because they provide the nodal structure of 
the trial function. With the fixed node approximation, the projected ground state has 
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the same nodal surfaces of the trial function, while the other details of the trial func- 
tion are automatically "optimized" for increasing projection time. It is thus important 
that the nodes provided by given spin-orbitals be accurate. Moreover, the optimiza- 
tion of nodal parameters (see below) is, in general, more difficult and unstable than 
for the pseudopotential parameters [6]. 

The simplest form of spin-orbitals for a system with translational invariance are 
plane waves (PW) 9k{r,a) = exp[i,k ■ r]. This form was used in the first QMC 
study of metallic hydrogen [33]. It is particularly appealing for its simplicity and 
still quahtatively correct since electron-electron and electron-proton correlations are 
considered through the "pseudopotential". The plane waves orbitals are expected to 
reasonably describe the nodal structure for metallic atomic hydrogen, but no infor- 
mation about the presence of protons appears in the nodes with PW orbitals. 

For insulating molecular hydrogen (i.e for Vg > 1-5), it is preferred to use lo- 
calized gaussian orbitals. There are different possibilities: a single isotropic gaussian 
centered at the middle of the bond was used in the first QMC study[33], while a 
single multivariate gaussian was used in the first CEIMC attempt[15, 16]. Another 
possibility is to form a molecular orbitals as Unear combination of two atomic gaus- 
sians orbitals centered on each proton: 

6k{r,Sk,i,Sk,2) = exp{-c\r - Sk,i\'^) + exp{-c\r - Sfe,2|^) (55) 

where Skj is the position of the j*'' proton of the A;*'' molecule. This kind of orbitals 
have a single variational parameter c. At present we are experimenting using a trial 
molecular wave function with these orbitals multiplied by the corrected RPA Jastrow 
within the CEIMC. 

The trouble with this strategy is that one should know which phase is stable be- 
fore performing the calculation. This is typical of ground state studies. However in 
CEIMC we would rather let the system find its own state for given temperature and 
density (or pressure). In particular, this approach is not appropriate to address the 
interesting region of molecular dissociation and metallization. This problem can be 
solved by using orbitals obtained as solution of a single-electron problem as in band 
structure calculations or in self-consistent mean field methods. In previous works on 
ground state hydrogen, the single electron orbitals for a given protonic state S, were 
obtained from a DFT-LDA calculation [34, 35, 36]. One of this study [35] estab- 
Ushed that energies from plane-waves determinants in metallic hydrogen are higher 
than the more accurate estimates from DFT-LDA orbitals by 0.05eV/atom at the 
density at which the transition between molecular and metallic hydrogen is expected 
(rg = 1.31). Obtaining the orbitals from a DFT-LDA calculation has, however, sev- 
eral drawbacks in coimection with the CEIMC. While for protons on a lattice we 
can solve the self-consistent theory for a primitive cell only, in a disordered con- 
figuration, we need to consider the entire simulation box. This is very expensive in 
CPU time and memory for large systems. Moreover, combining the LDA orbitals 
with Jastrow to improve the accuracy is not straightforward; substantial modification 
of the orbitals might be necessary requiring a reoptimization of the orbitals and the 
correlation factors, in principle, at each new ionic position. 
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Beyond the pair product trial action 

Over the years there have been important progress in finding trial functions substan- 
tially more accurate then the pair product form for homogeneous systems [12, 13]. 
Within the generalized Feynman-Kac formalism, it is possible to systematically im- 
prove a given trial function[13, 14]. The first corrections to the pair product action 
with plane wave orbitals are a three-body correlation term which modifies the corre- 
lation part of the trial function (Jastrow) and a "backflow" transformation which 
changes the orbitals and therefore the nodal structure (or the phase) of the trial 
function[14]. The new trial function has the form 

WTiR,S\S) = det[ek{x,,a,\S)]e-^'-^' (56) 

where U2 = J2i<j ^ij t^e two body "pseudopotential" discussed before, U3 the 
three-body term of the form 



t^3 = -E 



N 



(57) 



and finally the "quasiparticle" coordinates appearing in the plane wave orbitals are 
given by 

N 

Xi=ri + '^7]ij{rij)rij; (« = 1, • • • , iVe) (58) 

The functional form of the three-body term is that of a squared two-body force so its 
evaluation is not slower than a genuine 2-body term. In the homogeneous electron 
system, this term is particularly relevant at low density where correlation effects are 
dominant. On the other hand, the backflow transformation is more relevant at high 
density because in this limit the fermionic character of the system dominates. The 
same general framework should hold for hydrogen although a throughout study of 
the relative importance of those effects with density is still missing. The fundamental 
improvement of backflow orbitals for metallic hydrogen is that the nodal structure of 
the wave functions depends now on the proton positions. This provides better total 
energies for static protonic configurations [14] and improved energy differences and 
liquid structure in CEIMC [37, 16]. However it is not clear how appropriate this kind 
of wave function will be when entering in the molecular phase of hydrogen. 

The unknown functions, ^ij (r) , 77^ (r) in eqs. (57) and (58) need to be parame- 
terized in some way. In a first attempt we have chosen gaussians with variance and 
amplitude as new variational parameters [16]. This form was shown to be suitable for 
homogeneous electron gas [13]. Approximate analytical forms for (r) and i]ij{r), 
as well as for the two-body pseudopotential, have been obtained later in the frame- 
work of the Bohm-Pines collective coordinates approach [14]. This form is particu- 
larly suitable for the CEIMC because there are no parameters to be optimized. This 
trial function is faster than the pair product trial function with the LDA orbitals, has 
no problems when protons move around and its nodal structure has the same quality 
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as the corresponding one for the LDA Slater determinant [14]. We have extensively 
used this form of the trial wave function for CEIMC calculations of metalUc atomic 
hydrogen. 

Trial function optimization 

For metallic hydrogen we have described a parameter-free trial function which does 
not need optimization. However, if we use the pair proton action both for molecular 
or LDA orbitals, we are left with free parameters in the Jastrow factor and with the 
width of the gaussians for molecular orbitals. Optimization of the parameters in a 
trial function is crucial for the success of VMC. Bad upper bounds do not give much 
physical information. Good trial functions will be needed in the Projector Monte 
Carlo method. First, we must decide on what to optimize and then how to perform 
the optimization. There are several possibilities for the quantity to optimize and de- 
pending on the physical system, one or other of the criteria may be best. 

• The variational energy: Ey- If the object of the calculation is to find the least up- 
per bound one should minimize Ey. There is a general argument suggesting that 
the trial function with the lowest variational energy will maximize the efficiency 
of Projector Monte Carlo [38]. 

• The variance of the local energy: = J \H'1^\^ — Ey . If we assume that every 
step on a QMC calculation is statistically uncorrelated with the others, then the 
variance of the average energy will equal c^/p where p is the number of steps. 
The minimization of is statistically more robust than the variational energy 
because it is a positive definite quantity with zero as minimum value. One can 
also minimize a linear combination of the variance and the variational energy. 

• The overlap with the exact wave function: / If we maximize the overlap, 
we find the trial function closest to the exact wave function in the least squares 
sense. This is the preferred quantity to optimize if you want to calculate corre- 
lation functions, not just ground state energies since, then, the VMC correlation 
functions will be closest to the true correlation functions. Optimization of the 
overlap will involve a Projector Monte Carlo calculation to determine the change 
of the overlap with respect to the trial function so it is more complicated and 
rarely used. 

The most direct optimization method consists of running independent VMC calcula- 
tions using different set of numerical values for the variational parameters. One can 
fit the energies to a polynomial, performing more calculations near the predicted min- 
imum and iterating until convergence in parameter space is attained. The difficulty 
with this direct approach is that close to the minimum, the independent statistical 
errors will mask the variation with respect to the trial function parameters. This is 
because the derivative of the energy with respect to trial function parameters is very 
poorly calculated. Also, it is difficult to optimize, in this way, functions involving 
more than 3 variational parameters because so many independent runs are needed to 
cover the parameter space. 
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A correlated sampling method, known as reweighting [39, 40] is much more efficient. 
One samples a set of configurations {Rj} (usually several thousand points at least) 
according to some distribution function, usually taken to be the square of the wave- 
function for some initial trial function: \\I/t{R; {a}o)p. Then, the variational energy 
(or the variance) for trial function nearby in parameter space can be calculated by 
using the same set of points: 

Ev{a) = — — J , (59) 

where the weight factor, w{R) = {i'TiR', {«})/tf'T(^^; {"}o)|^, takes into account 
that the distribution function changes as the variational parameters change. One then 
can use a minimizer to find the lowest variational energy or variance as a function of 
{a} keeping the configurations fixed. However, there is an instability: if the param- 
eters move too far away, the weights span too large of a range and the error bars of 
the energy become large. The number of effective points of a weighted sum is: 

Neff = iY.w,r/Y.'vi (60) 

If this becomes much smaller than the number of points, one must resample and gen- 
erate some new points. When minimizing the variance, one can also simply neglect 
the weight factors. Using the reweighting method one can find the optimal value of 
wavefunction containing tens of parameters. 



2.5 Twist Average Boundary Conditions 

Almost all QMC calculations in periodic boundary conditions have assumed that the 
phase of the wave function returns to the same value if a particle goes around the pe- 
riodic boundaries and returns to its original position. However, with these boundary 
conditions, delocaUzed fermion systems converge slowly to the thermodynamic limit 
because of shell effects in the filling of single particle states. Indeed, with periodic 
boundary conditions the Fermi surface of a metal wiU be reduced to a discrete set 
of points in k-space. The number of k-points is equal to the number of electrons of 
same spin and therefore it is quite limited. 

One can allow particles to pick up a phase when they wrap around the periodic 
boundaries, 

^eiri + Ln, r2, • • •) = e'^^e{ri,r2, ■ ■ •)• (61) 

where we have assumed a cubic box of size L and n is a vector of integers. The 
boundary condition = is periodic boundary conditions (PBC), and the general 
condition with 6^0, twisted boundary conditions (TBC). If the periodic boundaries 
are used in all directions, each dimension can have an independent twist. Hence, in 
three dimension (3D), the twist angle is a three component vector. The free energy 
and therefore all equilibrium properties are periodic in the twist: F{d+2Trn) = F{0) 
so that each component of the twist can be restricted to be in the range — tt < 9i < tt. 
The use of twisted boundary conditions is commonplace for the solution of the band 
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structure problem for a periodic solid, particularly for metals. In order to calculate 
properties of an infinite periodic solid, properties must be averaged by integrating 
over the first Brillouin zone. 

For a degenerate Fermi liquid, finite-size shell effects are much reduced if the 
twist angle is averaged over: twist averaged boundary conditions (TABC)[30]. For 
any given property A the TABC is defined as 



TABC is particularly important in computing properties that are sensitive to the sin- 
gle particle energies such as the kinetic energy and the magnetic susceptibility. By 
reducing shell effects, accurate estimations of the thermodynamic Umit for these 
properties can be obtained already with a limited number of electrons. What makes 
this very important is that the most accurate quantum methods have computational 
demands which increase rapidly with the number of fermions. Examples of such 
methods are exact diagonalization (exponential increase in CPU time with N), vari- 
ational Monte Carlo (VMC) with wave functions having backflow and three-body 
terms [13] (increases as iV^), and transient-estimate and released-node Diffusion 
Monte Carlo methods [41] (exponential increase with N). Moreover, size extrapo- 
lation is impractical within CEIMC since it would have to be performed for any 
proposed ionic move prior to the acceptance test. Methods which can extrapolate 
more rapidly to the thermodynamic limit are crucial in obtaining high accuracy. 

Twist averaging is especially advantageous in combination with stochastic meth- 
ods (i.e. QMC) because the twist averaging does not necessarily slow down the 
evaluation of averages, except for the necessity of doing complex rather than real 
arithmetic. In a metallic system, such as hydrogen at very high pressure, results in 
the thermodynamic limit require careful integration near the Fermi surface because 
the occupation of states becomes discontinuous. Within LDA this requires "k-point" 
integration, which slows down the calculation linearly in the number of k-points re- 
quired. Within QMC such k-point integration takes the form of an average over the 
(phase) twist of the boundary condition and can be done in parallel with the aver- 
age over electronic configurations without significantly adding to the computational 



In CEIMC we can take advantage of the twist averaging to reduce the noise in 
the energy difference for the acceptance test of the penalty method (see below). In 
the electron gas, typically 1000 different twist angles are required to achieve conver- 
gence [30]. We have used the same number of twist angles in CEIMC calculations of 
metallic hydrogen. Different strategies can be used to implement the TABC [30]. We 
have used a fixed 3D grid in the twist angle space, at each grid point run independent 
QMC calculations and then averaged the resulting properties. This procedure can be 
easily and efficiently implemented on a parallel computer. Recently, we have devised 
a sampUng procedure to randomize the grid points at each ionic step. We have lim- 
ited experience, but, so far, we have evidence that good convergence in the electronic 
and ionic properties can already be reached for a number of twist angles as low as 30 




(62) 



effort. 



[42]. 
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2.6 Sampling electronic states: the "bounce" algorithm 

In this section we describe the way we have implemented electronic move in the 
CEIMC method. In particular we present an original algorithm for RQMC, particu- 
larly suitable for CEIMC, called the "bounce" algorithm. 

First, how do the particles move in VMC? In the continuum it is usually more 
efficient to move the particles one at a time by adding a random vector to a parti- 
cle's coordinate, where the vector is either uniform inside of a cube, or is a normally 
distributed random vector centered around the old position. Unfortunately, this pro- 
cedure cannot be used with backflow orbitals. This is because the backflow trans- 
formation couples all the electronic coordinates in the orbitals so that once a single 
electron move is attempted the entire Slater determinant needs to be recomputed, an 
0{N^) operation. It is much more efficient to move all electrons at once, although 
global moves could become inefficient for large systems. 

Next we describe the electronic sampling within RQMC. In the original work on 
RQMC [18], the electronic path space was sampled by a simple reptation algorithm, 
an algorithm introduced to sample the configurational space of linear polymer chains 
[43]. Remember that in RQMC the electronic configurational space is the space of 
SA^e-dimensional random paths of length "t". In practice, the imaginary time is dis- 
cretized in M time slices t = t/M and the paths become discrete linear chains 
of M + 1 beads. Let us indicate with Q the entire set of 3N{M + 1) coordinates 
Q = {Ro, • • • , Rm}- According to eqs. (33), (34) and (36), the path distribution is 

71(g) = |iZ'r(i?o)>Z'T(i?M)|e"-^^=i^'^'"'=-^'^'='"^e"^L^^+^^=i ^=^(-R'=)+^^J 

(63) 

Given a path configuration Q, a move is done in two stages. First one of the two ends 
(either Rq or Rm) is sampled with probability 1/2 to be the growth end Rg. Then a 
new point near the growth end is sampled from a Gaussian distribution with center 
at Rg + 2\TeF(Rg). In order to keep the number of links on the path constant, the 
old tail position is discarded in the trial move. The move is accepted or rejected with 
the MetropoUs formula based on the probability of a reverse move. For use in the 
following, let us define the direction variable d as d = +1 for a head move (Rg = 
Rm), and d = —1 for a tail move (Rg = Rq). In standard reptation, the direction d 
is chosen randomly at each attempted step. The transition probabiUty P{Q — > Q') 
is the product of an attempt probability Td{Q Q') and an acceptance probability 
0'd{Q Q')- The paths distribution n{Q) does not depend on the direction d in 
which it was constructed. In the Metropolis algorithm, the acceptance probability for 
the attempted move is 



ad{Q ^ Q') 



^ n{Q')T^,{Q' ^ Q) 



(64) 



n{Q)TdiQ ^ Q') 

which ensures that the transition probability Pd{Q — > Q') satisfies detailed balance 

n{Q)Pd{Q ^ Q') = n{Q')P_d{Q' ^ Q) (65) 
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The autocorrelation time of this algorithm, that is the number of MC steps be- 
tween two uncorrelated configurations, scales as [M'^ /A], where A is the acceptance 
rate of path moves, an unfavorable scaling for large M (i.e. large projection time 
t). Moreover the occasional appearance of persistent configurations bouncing back 
and forth without really sampling the configuration space has been previously ob- 
served [44]. These are two very unfavorable features, particularly in CEIMC, where 
we need to perform many different electronic calculations. There is a premium for a 
reliable, efficient and robust algorithm. 

We have found that a minimal modification of the reptation algorithm solves 
both of these problems. The idea is to chose randomly the growth direction at the 
beginning of the Markov chain, and reverse the direction upon rejection only, the 
"bounce" algorithm. 

What follows is the proof that the bounce algorithm samples the correct prob- 
ability distribution 11 (Q). The variable d is no longer randomly sampled, but, 
as before, the appropriate move is sampled from the same Gaussian distribution 
Td{Q — > Q') and accepted according to the Eq. (64). In order to use the techniques 
of Markov chains, we need to enlarge the state space with the direction variable d. 
In the enlarged configuration space {Q, d}, let us define the transition probabihty 
P{Q, d — > Q' , d') of the Markov chain. The algorithm is a Markov process in the 
extended path space, and assuming it is ergodic, it must converge to a unique station- 
ary state, T(<5, d) satisfying the eigenvalue equation: 

^ r(g, d) P{Q, d ^ Q', d') = T{Q', d'). (66) 

Q,d 

We show that our desired probability n{Q) is solution of this equation. Within the 
imposed rule, not all transitions are allowed, but P{Q, d Q\ d') ^ for d = d' 
and Q ^ Q' (accepted move), or d! = —d and Q = Q' (rejected move) only. Without 
loss of generality let us assume d' = +1 since we have symmetry between ±1. Eq. 
(66) with Y{Q, d) replaced by iT(Q) is 

n{Q')P{Q', -1 ^ Q', i)+Yl n{Q)P{Q, i - Q', i) = n{Q'). (67) 

Because of detailed balance Eq.(65), we have 

niQ)P{Q, 1 ^ Q', 1) = n{Q')P{Q', -1 ^ Q, -1) 
which, when substituted in this equation gives 



n{Q') 



P{Q', -1 ^ Q', 1) + ^ P(Q', -1 ^ Q, -1) 
Q 



n{Q'). (68) 



Note that we have completed the sum over Q with the term Q = Q' because its 
probabihty vanishes. The term in the bracket exhausts all possibihties for a move 
from the state (Q', —1), thus it adds to one. Hence n{Q) is a solution of eq. (66) 
and by the theory of Markov chains, it is the probability distribution of the stationary 
state. 
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3 Sampling ionic states: The penalty method 

In Metropolis Monte Carlo a Markov chain of ionic states S is generated accord- 
ing to the Boltzmann distribution P{S) oc e'^^'""^^^ where Ebo{S) is the Bom- 
Oppenheimer energy for the ionic configuration S and f3 the inverse temperature. 
From the state S a trial state S' is proposed with probabihty T{S S') and the 
detailed balance condition is imposed by accepting the move with probabihty 



A{S ^ S') = min 



T{S' ^ S)e-^^Bo(s') 
' T{S ^ S')e-^^Bo(S) 



(69) 



Under quite general conditions on the system and on the a-priori transition proba- 
bihty T{S S'), after a finite number of MC steps the Markov chain so generated 
will visit the states of the configurational space with a frequency proportional to their 
Boltzmann's weight [43]. 

In CEIMC estimate of Ebo (S) is affected by statistical noise. If we ignore the pres- 
ence of noise and we use the standard Metropolis algorithm, the results will be bi- 
ased, the amount of bias increasing for increasing noise level. A possible solution 
would be to run very long QMC calculations in order to get a negligibly small noise 
level resulting in a negligible bias. However, the noise level decreases as the number 
of independent samples to the power 1/2, that is to decrease the noise level by one 
order of magnitude we should run 100 times longer, an unfavorable scahng if we re- 
alize that we have to repeat such calculation for any attempted move of the ions. The 
less obvious but far more efficient solution is to generalize the Metropolis algorithm 
to noisy energies. This is done by the Penalty Method[45]. The idea is to require the 
detailed balance to hold on average and not for any single energy calculation. 

Let us consider two ionic states {S, S') and call d{S, S') the "instantaneous" en- 
ergy difference times the inverse temperature. Let us further assume that the average 
and the variance of 6{S, S') over the noise distribution P{6\S S') exhist 

A{S,S') = P]Ebo{S') - Ebo{S)] =< 5{S,S') > (70) 
(7^(5, S') = < {S{S, S') - A{S, S')f > (71) 

We introduce the "instantaneous" acceptance probabihty, a{6\S, S') and impose the 
detailed balance to hold for the average of a(^|S', S') over the distribution of the 
noise 

A{S S') = e-^^^'^'^A{S' S) (72) 

where 

/oo 
ddP{6\S ^ S')a{S\S,S') (73) 
-oo 

and we have defined 

'T{s' ^sy 



r{S, S') = A{S, S') - In 



T{S- 



(74) 
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If we assume the quite general conditions a{S\S, S') = a{6) and P{6\S S') = 
P{—6\S' —^ S) to hold, the detailed balance can be written 

/oo 
d6P{6\S S') [a{S) - e-^a{-S)] = (75) 
-oo 

which, supplemented by the condition a{S) > 0, is the equation to solve in order to 
obtain the acceptance probabiUty a{S). The difficulty is that during the MC calcula- 
tion we do not know either P{S\S S') nor A{S, S'). 

In order to make progress let us assume, as it happens in many interesting cases, that 
the noise of the energy difference is normally distributed so that 

P{d\S ^ S') = (27r£72)-V2 exp [-{S - Af/2a^] (76) 

Let us, moreover, assume that we know the value of the variance a. It is not difficult 
to check that the solution of eq. (75) is 



an{S\<j) = min 



T{S' ^ S) „ 



(77) 



The uncertainty in the energy difference just causes a reduction in the acceptance 
probability by an amount exp(— (t2/2) for 5 > —a^/2. The integral of a„(5|cr) over 
the guassian measure provides 

A{S ^ 5') = i crfc { [ay 2 + r{S ^ S')] /(V2(t)} + 

i erfc { [aV2 - r{S ^ S')] /{V2a) } e'^^^^^'^ (78) 

which satisfies eq. (72) since r{S' ^ S) = -r(5 ^ S'). Note that eq. (77) 
reduces to the standard Metropolis form for vanishing a [43]. 

An important issue is to verify that the energy differences are normally dis- 
tributed. Recall that if the moments of the energy difference are bounded, the central 
limit theorem implies that given enough samples, the distribution of the mean value 
will be Gaussian. Careful attention to the trial function to ensure that the local ener- 
gies are well behaved may be needed. 

In practice not only the energy difference A but also the variance a is unknown 
and must be estimated from the data. Let us assume that, for a given pair of ionic 
states {S, S'), we generate n statistically uncorrelated estimates of the energy dif- 
ference {yi. . . . , y,,} each normally distributed with their first and second moments 
defined in the usual way 

A=<yi> (79) 
a^ = < ivi - Af > (80) 



Unbiased estimates of A and are 
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1 " 

^ = - V yi 



n . 

1=1 



(81) 



n(n-l)^ 

An extension of the derivation above provides [45] 

T(S" ^ 5) 



a{5, X ,n) = min 



T(5- 



exp(— 5 — ub) 



(83) 



where 



^^ = T + 4(;^+ 3(n+l')(n + 3) +--- ^'^^ 

The first term in eq. (84) is the penalty in the case we know the variance. The error 
of the noise causes extra penalty which decreases as the number of independent sam- 
ples n grows. In the limit of large n the first term dominates and we recover eq. (77). 
Eq. (84) must be supplemented by the condition x^/'^ < 1/4 for the asymptotic ex- 
pansion (84) to converge and the instantaneous acceptance probabiUty to be positive 
[45]. 

The noise level of a system can be characterized by the relative noise parameter, 
/ = (/3cr)^i/t(), where t is the computer time spent reducing the noise, and to is the 
computer time spent on other pursuits, such as optimizing the VMC wave function or 
equilibrating the RQMC runs. A small / means little time is being spent on reducing 
noise, where a large / means much time is being spent reducing noise. For a double 
well potential, the noise level that gives the maximum efficiency is around Pa w 
1, with the optimal noise level increasing as the relative noise parameter increases 
[45]. In CEIMC runs for hydrogen the noise level pa ranges between 0.3 and 3, the 
optimal value being around 1. 



3.1 Efficient strategies for electronic energy differences 

As explained above, we need to evaluate the energy difference and the noise between 
two protonic configurations (5, S'). The distance in configurational space between 
S and S' in an attempted move is however quite limited since we have to move 
all protons at once. Indeed, each backflow orbital depends on the position of all 
protons and single proton moves would require recomputing the entire determinant 
for each attempted move, a 0(7V^) operations for a global move. Instead, moving all 
protons together requires a single determinant calculation per ionic move, a O(iV^) 
operation. In this case, performing two independent electronic calculations for 5 
and S' to estimate the energy difference would be very inefficient since AEbo « 
Ebo{S) and the noise on the energy difference would just be twice the noise on 
the single energy estimate. The strategy to adopt is to compute the energy difference 
from correlated sampling, i.e. sampling the electronic configurational space from a 
distribution which depends both on S and S' and estimating the energy difference 
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and the other electronic properties of interest by reweighting the averages [16]. It is 
possible to show that the optimal sampling function, i.e. the sampling distribution 
for which the variance of the energy difference is minimal, takes the form 

P{Q\S,S') oc \n{Q\S){EBo{S)- < Ebo{S) >) - n{Q\S'){EBo{S')- < Ebo{S') >)| 

(85) 

where < Ebo > is the estimate of the BO energy. In order to use this sampUng 
probability we need to estimate the BO energies of the two states before performing 
the sampling. A simpler form which avoid this problem is 

P{Q\S,S')(xn{Q\S) + n{Q\S') (86) 

We emphasize that the "reptile" space for the electron paths depends on the proton 
coordinates so that, because of the fixed-node restriction for fermions, legal paths for 
S may or may not be legal for S' and vice-versa. These two forms of importance 
sampling have the property that they sample regions of both configuration spaces (S 
and S') and make the energy difference rigorously correct with a bounded variance. 



3.2 Pre-rejection 

We can use multi-level sampling to make CEIMC more efficient [19]. An empirical 
potential is used to "pre-reject" moves that would cause particles to overlap and 
be rejected anyway. A trial move is proposed and accepted or rejected based on a 
classical potential 



Ai = min 



(87) 



where AVd = Vci{S') — Vci{S) and T is the sampling probability for a move. If it 
is accepted at this first level, the QMC energy difference is computed and accepted 
with probability 

A2 = min [1, exp(-/3Z\£;BO - «s) exp(/3Z\V;0] (88) 

where us is the noise penalty. Compared to the cost of evaluating the QMC energy 
difference, computing the classical energy difference is much less expensive. Reduc- 
ing the number of QMC energy difference evaluations reduces the overall computer 
time required. 

For metallic hydrogen a single yukawa potential was always found to be suit- 
able for pre-rejection. For molecular hydrogen we use instead a Silvera-Goldman 
potential [46] riparametrized in such a way to have the center of interaction on any 
single proton in the molecule rather than on the molecular center of mass. In prac- 
tice, we take a single Yukawa potential for intermolecular proton-proton interaction 
and the Kolos-Wolniewski potential [47] for the bonding interaction. The Yukawa 
screening length and the prefactor are optimized to reproduce the results of the 
Silvera-Goldman model. This new potential is suitable for pre-rejecting all types 
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of moves that we attempt in the molecular hydrogen, namely molecular rotations, 
bond stretching and molecular translations. The original Silvera-Goldmann potential 
being spherically symmetric around the molecular center is not suitable to pre-reject 
the rotational moves. 



4 Quantum protons 

By increasing pressure and/or decreasing temperature, ionic quantum effects can be- 
come relevant. Those effects are important for hydrogen at high pressure [7, 48 J. 
Static properties of quantum systems at finite temperature can be obtained with the 
Path Integral Monte Carlo method (PIMC) [19]. We need to consider the ionic ther- 
mal density matrix rather than the classical Boltzmann distribution: 

Pp{S,S'\/3) =< S\e-l^^^''+^^'^^\S' > (89) 

where Kp is the ionic kinetic energy operator and /3 = {KbT)~^ is the inverse 
physical temperature. Thermal averages of ionic operators Ap (diagonal in configu- 
rational space) are obtained as 

M0) = Y{0)j '^^MS)Pvi.s, s\0) (90) 

where Z is the partition function 

Z{I3) = j dSpp{S,S\l3)=e-^^ (91) 

and F is the Helmholtz free energy of the system. As before, the thermal density 
matrix can be computed by a factorization of /3 in many (?) small intervals (time 
slices Tp = (3/P) and by a suitable approximation for the "high temperature" (or 
short time) density matrix. According to the Feynman-Kac formula (see eqs. (24) 
and (35)) the diagonal part of the thermal density matrix is the sum over all closed 
paths, i.e. paths starting at S and returning to S after a "time" /?. This is the famous 
"isomorphism" between quantum particles and ring polymers [21, 19]. Considering 
particle statistics in the PIMC is more difficult than in RQMC. The reason is that in 
PIMC the state of the classical system, which has no symmetry built in, plays the 
role of the trial functions in RQMC and permutations need to be sampled explic- 
itly with an additional level of difficulty in the method [19, 26]. However, protonic 
statistics become relevant when the quantum dispersion is comparable to the interi- 
onic distances Ap = y/2Xpl3 « {Np/V)~^/^ = Up^^^. This define a degeneracy 

temperature kBToinp) = 2\pr^J^ below which quantum statistics need to be con- 
sidered. For hydrogen Td — 66.2(X)/r^, where Vg is the usual ion sphere radius of 
coulomb systems Vg = (3/47rnj,)^/^. Therefore proton statistics in metallic hydro- 
gen {Vs < 1.3) becomes relevant below 50K depending on the density, a regime that 
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we have not investigated yet. Proton statistics in molecular hydrogen is also quite im- 
portant and results in the separation between ortho- and para-hydrogen[49]. Because 
this effect is relevant only at low temperature, we have disregarded it as well. 

In order to implement the PIMC we need a suitable approximation for the high 
temperature density matrix Pp{S, S'\ti,). We could use either the primitive approxi- 
mation or the importance sampling approximation described earher. However a bet- 
ter approximation, in particular for distinguishable particles, is the pair product ac- 
tion [19] which closely resembles the pair trial function. The idea is to build the 
many body density matrix as the product over all distinct pairs of a two-body density 
matrices obtained numerically for a pair of isolated particles. At high temperature the 
system approaches the classical Boltzmann distribution which is indeed of the pair 
product form. The method is described in detail in ref. [19]. Here we just explain 
how we can take advantage of this methodology within the CEIMC scheme. In order 
to use the method of pair action, we need to have a pair potential between quan- 
tum particles. In CEIMC, however, the interaction among protons is provided by the 
many-body BO energy. Our strategy is to introduce an effective two-body potential 
between protons Ve and to recast the ionic density matrix as 

p^{S,S'\Tp) =< S\e-^^^"^+^^''°-^'^^S' > 

w< S\e-''''"'\S' > e-^[^«o(S)-Ve(S)]+[BBo(s')-Ve(s')] (92) 

where He = Kp + Ve and the corrections from the effective potential to the true 
BO energy are treated at the level of the primitive approximation. We can compute 
numerically the matrix elements of the effective pair density matrix 0\tp) as ex- 
plained in ref. [19]. The effective A'p-body density matrix is approximated by 

Np 
ij 

(93) 

where po is the free particle density matrix and Ue{sij, •s^^ Itj,) is the effective pair 
action. The explicit form for the partition function is then 

Z{P) = [dS,...dSp f[ exp I -^-^^-^^ - ^«e(4> 
fc=i [ ^'^P^P ij 

exp j-Tp YjiEBoiSk) - VeiSk)]^ (94) 

with the boundary condition: Sp+i = 5i. As for the pre-rejection step in the pro- 
ton moves, we have used different effective potentials according to the system un- 
der consideration. In metalhc hydrogen (a plasma) we have used a smooth screened 
coulomb form and found that it provides a fast convergence in Tp. Convergence can 
be assessed by monitoring the various terms in the estimator for the proton kinetic 
energy[19]. A good effective potential should provide uniform convergence of the 
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various orders in Tp. With this effective potential, we have found convergence to the 
continuum limit (jp 0) for 1/rp > SOOOX which allows to simulate systems at 
room temperature with only M 10 proton slices for > 1. In molecular hy- 
drogen we need to consider the extra contribution of the bonding potential. We have 
used the Kolos-Wolniewicz bonding potential in connection with the same smooth 
screened coulomb potential for the non bonding interactions. Convergence with Tp is 
observed in a similar range. 

A very nice feature of ionic PIMC in CEIMC is that considering ionic paths rather 
then classical point particles does not add any computational cost to the method. Let 
us suppose we run classical ions with a given level of noise {f3aci)'^. Consider now 
representing the ions by P time slices. To have a comparable extra-rejection due to 
the noise we need a noise level per shce given by: (rpCTft)^ ~ {fiaciY /P which 
provides ~ Pcr^i . We can allow a noise per time sUce P times larger which means 
considering P times less independent estimates of the energy difference per slice. 
However we need to run P different calculations, one for each different time shce, so 
that the amount of computing for a fixed global noise level is the same as for classical 
ions. In practice, however, because our orbitals depend on all proton positions, we 
are forced to move the proton positions at given imaginary time all together with a 
local (in imaginary time) update scheme. It is well known that the autocorrelation 
time of schemes with local updates rapidly increases with the chain length and this 
is the ultimate bottleneck of our present algorithm [19]. It is therefore essential to 
adopt the best factorization in order to minimize the number of time slices P needed 
and therefore the efficiency of the method. 

When using TABC with quantum ions, for any proton time shce we should, in 
principle, perform a separate evaluation of the BO energy difference averaged over 
all twist angles. Instead at each protonic step, we randomly assign a subset of twist 
phases at each time slice and we compute the energy difference for that phase only. 
The TABC is then performed by adding up aU the contributions from the different 
time slices. We have checked in few cases that this simpUfied procedure does not 
give detectable biases in the averages. 

In practice we move all sUces of all protons at the same time by a simple random 
move in a box and we pre-reject the moves with the effective pair action. 

5 Summary of results on high pressure hydrogen 

In this section we briefly summarize some of the CEIMC results we have obtained for 
high pressure hydrogen. Figure 1 shows what is known and what has been predicted 
about the hydrogen phase diagram in a wide range of pressures and temperatures. The 
rectangular region in the right upper comer is the region where R-PIMC method can 
make reliable predictions [50, 51, 52]. There have been many studies of the ground 
state of hydrogen (T=0) including some QMC investigations [33, 34, 35, 36]. They 
have predicted a metalhzation density corresponding at = 1.31 accompanied by 
a transition from a molecular m-hcp structure to a diamond lattice of protons. At 
intermediate temperatures a number of ab-initio Molecular Dynamics studies have 
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Fig. 1. Hydrogen phase diagram. Continuous transition lines are experimental results, dashed 
lines are theoretical prediction from various methods. Squares and right-triangle are ab-initio 
MD predictions of molecular melting [57] and molecular dissociation in the liquid phase [56]. 
The diamonds are shock-waves experimental data through the liquid metalization [58]. The 
triangles are earlier CEIMC data for the insulating molecular state [15, 16] while the green 
domain on the extreme right indicates the CEIMC prediction for the melting [17]. Red lines 
are model adiabats for the interior of the giant planets of the solar system [59]. 

been performed in the molecular and in the metallic phases, both in the crystal and 
in the liquid state [53, 54, 55, 56, 57]. 

Metallic hydrogen 

We first focus on the metallic system for pressure beyond the molecular dissociation 
threshold. In this region, hydrogen is a system of protons and delocalized electrons. 
At low enough temperature the protons order in a crystalline lattice which melts upon 
increasing temperature. The low temperature stable structure as a function of density 
is still under debate. The most accurate ground state QMC calculation [35], indicates 
that hydrogen at the edge of molecular dissociation will order in a diamond struc- 
ture, and upon increasing density will undergo various structural transformations 
ultimately transforming to the bcc structure. However, these prediction are extrapo- 
lated from a single calculation at — 1.31 and temperature effects are absent. With 
CEIMC we have investigated the density range G [0.8, 1.2] and the temperature 
range T G [SOO/C, SOOOiC] across the melting transition of the proton lattice and 
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up to the lower limit of applicability of RPIMC. We limited the study to systems of 
Np — 32 and Np — 54 which, for cubic simulation boxes, form fee and bcc lattices, 
respectively. We have observed the melting and refreezing of the protons and made 





2 2 

a I H atom 



Fig. 2. Ne ~ Np — 16, Vs — 1.31. Dependence of total energy, variance and energy differ- 
ence for a pair of proton configurations {S, S') on the RQMC projection time. The study is 
performed for Te = 0.02H~^ . Dotted lines represent the variational estimates with their error 
bars. In panel b) and c) the lines are exponential fits to data and in panel d) the continuous 
line is a linear fit in the region < 0.005. Black circles (3BF-A) are results obtained with 
the analitical three-body and backflow trial wave functions discussed earlier, the red triangle 
is a variational result with a Slater- Jastrow trial function with simple plane wave orbitals and 
the blue squares are results from a trial function with LDA orbitals and an optmized two-body 
Jastrow. 

A number of interesting question about the convergence and the efficiency of 
the CEIMC algorithm need to be answered before starting a systematic study. An 
important one is: how large must the electronic projection time be in order to get 
convergence in the energy difference to the ground state value and therefore obtain 
unbiased sampling in CEIMC? In order to answer such a question we have selected 
a pair of protonic configuration {S, S') at given density and computed the energy 
difference, together with total energy and variance, versus the electronic projection 
time. The results reported in the figure 2 correspond to a system of = Np = 16 
at rs — 1.31 with the twist phase 6 = (0.4, 0.5, 0.6)7r. Results are obtained for an 
electronic time step Te = 0.02H^^, a compromise between accuracy and efficiency. 
In panel a) we report the energy difference versus the projection time t to show 
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that it does not depend on the projection time when using the accurate trial func- 
tion with 3-body terms and backflow orbitals (3BF-A, black dots) discussed above. 
This suggests that the proton configuration space can be sampled using VMC for the 
electrons which is faster and more stable than RQMC. On the same panel we have 
reported the VMC estimate obtained with a trial function with a Slater determinant 
of simple plane waves and a two-body RPA Jastrow (SJ, red triangle) and a RQMC 
result obtained for a trial function with a 2-body RPA-Jastrow and a Slater determi- 
nant of self consistent LDA orbitals (LDA-J, blue squares). The simple SJ function 
at the variational level has a much larger energy difference (in absolute value) and 
therefore will provide a biased sampling of protonic configurations. It is not clear 
whether the RQMC projection with such trial function will recover the correct en- 
ergy difference (remember that the two kind of trial functions have different nodal 
structure). On the other hand, no difference is detected between the RQMC results 
from the 3BF-A and the LDA-J trial functions at the same projection time. This 
gi • - ■ - . 
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Fig. 3. Ne ^ Np = 54, = 1.2, T = SOOOir, T point. Comparison of VMC (black) 
and RQMC (red) data for the pair correlation functions. The projection time in RQMC is 
te = 0.68H~^ with an electronic time step of Te = O.OIH^^. Protons are considered as 
classical point particles. 

the variance of configuration S versus t are reported in panels b) and c), respec- 
tively. The two panels nicely illustrate the limiting process toward the fixed-node 
ground state operated by the projection (see eq.(l 1) ). Within the 3BF-A trial func- 
tion we observe a large gain in energy with the projection time, the difference being 
£;(oo) - E{0) = 5.7mH/at = 1810K/at. In both panels the results of the SJ tiial 
function are off scale, while we have reported the result of the LDA-J wave function. 
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According to the total energy, the quality of LDA-J function is superior to the one 
of 3BF-A function for this configuration. The same quality was instead detected for 
the bcc lattice configuration [14]. Finally, panel d) we report the total energy versus 
the variance of the 3BF-A trial function to show that both quantities go linearly with 
the error in the wave function when we are close enough to the ground state (the 




Fig. 4. Ne ^ Np = 54, = 1.0, T = lOOOA', T point. Comparison of VMC (black) and 
RQMC (red) data for the pair correlation functions. The projection time in RQMC is — 
l.OH^^ with an electronic time step of Te — 0.02_ff^^. Protons are considered as classical 
point particles. The proton-proton pair correlation functions exhibits a structure reminiscent 
of the bcc crystal structure. RQMC data for Qepir) and gl^{r) show a finite time step errors 
at short distances which however do not contribute significantly to the energy. 

with VMC energy differences, we have performed two test runs for the systems of 
Np ^ Ne ^ 54 at the F point (zero twist angle), one at = 1.2 and T = 5000 A' 
and the other at = 1 and T = lOOOX. The comparison for the pair correlation 
functions are shown in figures 3 and 4 respectively. 

It is interesting to compare the predictions of CEIMC with other methods. In 
refs. [16, 17] we have compared with Restricted Path Integral Monte Carlo data and 
with Car-Parrinello Molecular Dynamics data with LDA forces (CPMD-LDA) [54]. 
In figure 5 we compare CEIMC and CPMD-LDA 5pp(r)'s for classical protons at 

— 1 and various temperatures. Both calculations are done with PBC (zero twist 
angle) and CEIMC uses 54 protons while CPMD-LDA used 162 protons (both are 
closed shells in the reciprocal space, so that the electronic ground state is not de- 
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Fig. 5. Proton-proton pair correlation functions for classical protons at = 1 and various 
temperatures. Comparison between CEIMC and CPMD-LDA data. 



generate). We see that in a wide range of temperatures, CPMD-LDA and CEIMC 
are off by a factor of 2 in temperature. CPMD-LDA predicts a less structured fluid 
and locates the melting transition at roughly T„j ~ 350K [54]. With CEIMC instead 
more structure is found and the melting is located between lOOOK and 1500K [17]. 
Though several reasons could be at the origin of such unexpected discrepancy, we 
believe that the probelm arises from a too smooth BO energy surface as provided 
by LDA. Evidences of this fact are also provided by previous ground state calcula- 
tions for hydrogen crystal structures [35] which similarly found energy differences 
between various structures from LDA to be roughly half of the coiTesponding one 
from DMC. This will explaine the factor of 2 in the temperature scale observed in 
figure 5. 

In ref. [17] we have reported data for the equation of state of metallic hydrogen. 
Proton quantum effects are quite important at such high density and need to be con- 
sidered carefully to get accurate prediction of the equation of state. The importance 
of proton quantum effect is partially reported in ref. [60]. Also we have shown that 
electronic VMC with 3BF-A trial function is accurate enough to sample the proton 
configuration space. However RQMC should be used in order to obtain accurate re- 
sults for the energy and the pressure. A good strategy is to run CEIMC with VMC to 
sample efficiently the proton configurations and then run RQMC on fixed statistical 
independent protonic configurations previously generated. 
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Insulating Molecular Hydrogen 

In this section we report very preliminary results in the insulating molecular liquid 
phase. We have investigated the state point at = 2.1 and T = 4530K because 
at this point, and other scattered points around it, experimental data for the equation 
( 
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Fig. 6. Comparison between VMC and RQMC computed pair correlation functions for insu- 
lating molecular hydrogen at = 2.1, T = 45307^. Np = Ne = 54. 

[15, 16] focus on the same point although using quite a different trial functions. We 
used a pair product trial function with modified RPA Jastrows (electron-electron and 
electron-proton) and Slater determinants built with the molecular orbitals of eq. (55). 
At variance with the metallic trial function, we have now 5 variational parameters to 
be optimized, 4 in the Jastrow factor and a single one in the orbitals. In particu- 
lar, the nodal structure of the wave function will be affected by the width of the 
gaussian orbitals used to build the molecular orbitals. In the first CEIMC attempt 
[15, 16], a single multivariate gaussian, centered at the molecular center of mass and 
with a width different for each molecule was used. Thus the number of variational 
parameter was equal to three time the number of molecules. Also, the entire opti- 
mization procedure for each attempted protonic configuration was performed prior 
the Metropolis test. Clearly, the optimization step was a bottleneck of that scheme. 
In our present scheme, we have gained evidence that we can optimize the parameters 
on a single protonic configuration (even the parameters optimized on a lattice config- 
uration are good) and use the same values for simulating the liquid phase. A similar 
conclusion was obtained in the metallic phase with numerically optimized orbitals 
[16]. In this way the optimization step needs to be performed only upon changing 
density and/or number of particles, while we can use the same set of values for the 
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variational parameters to span the temperature axis at fixed density and number of 
particles. In figure 6 we compare the pair correlation functions as obtained by VMC 
and RQMC with = O.Olif-^ and = lAH'^ for a system of 27 hydrogen 
molecules. Since the system is insulating we do not average over the twist angle, 
but just use the F point. We observe that RQMC enhances slightly the stregth of 
the molecular bond and the electronic correlation inside molecules with respect to 
VMC. Good agreement with the early CEIMC results is observed not only for the 
correlation function, but also for the equation of state. In our present calculation we 
obtain P[RQMC) = 0.224(5)M6ars to be compared with the previous estimate 
of 0.225(3)M6ars [15, 16] and with the experimental data P = Q.2MMbars. The 
deviation from the experimental data is only roughly 5% and it is particularly en- 
couraging if we consider that our data are for classical protons and quantum effects 
are expected to sUghtly increase the pressure. 

6 Conclusions and future developments 

In this paper we have described the principles of CEIMC and given some technical 
details on its practical implementation. The results for metallic and molecular hydro- 
gen show that CEIMC is a practical strategy to couple ground state QMC methods 
for the electronic degrees of freedom with a finite temperature Monte Carlo simula- 
tion of the ionic degrees of freedom. To our knowledge, CEIMC is the only method 
available so far to perform ab-initio simulations based on QMC energies. 

In a recent work, Grossman and Mitas have proposed a strategy to correct ab- 
initio Molecular Dynamics energies with QMC [62]. This is however different from 
CEIMC because the nuclear degrees of freedom are still sampled on the basis of DFT 
forces. Therefore, when applied to metalhc hydrogen, that method would have found 
the same Uquid structure as obtained by CPMD-LDA [54]. 

Very recently an interesting proposition on how to use noisy forces in ab-initio 
Molecular Dynamics has appeared [63]. The same strategy could be used with noisy 
QMC forces to simulate the dynamics of classical ions. A first attempt has already 
shown the feasibilty of this promising method [64] although the results are stiU very 
preliminary. 

A crucial aspect of CEIMC is the choice of the electronic trial wave function. 
The ones we have discussed in this paper are either suitable for the metallic state 
or for the molecular state and their quahty in describing the metalization-molecular 
dissociation transition is questionable. A current development within CEIMC is an 
efficient strategy to generate on the fly single electron orbitals depending on the 
instantaneous ionic configuration, in the spirit of the LDA orbitals previously used 
[35] . We have devised an efficient algorithm which is able to provide accurate orbitals 
in a reasonable computer time. We are at present using these orbitals to explore the 
molecular dissociation under pressure in the liquid state [65]. The same strategy 
can be used to obtain accurate prediction for the hydrogen equation of state in the 
temperature range between 50K and 5000K. Also application to helium and heUum- 
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hydrogen mixtures, very interesting systems for planetary physics, can be envisaged 
with the same methodology. 

Tests for non-hydrogenic systems are needed to find the performance of the 
method on a broader spectrum of apphcations. The use of pseudopotentials within 
QMC to treat atoms with inner core is well tested [6]. What is not clear is how much 
time will be needed to generate trial functions, and to reduce the noise level to accept- 
able hmits. Clearly, further work is needed to allow this next step in the development 
of microscopic simulation algorithms. 
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